SOCP 


Software  for 

Second-order  Cone  Programming 

User’s  Guide 


Beta  Version  April  1997 


Miguel  Sousa  Lobo 
Lieven  Vandenberghe 
Stephen  Boyd 


mloboOisl . stanford.edu 
vandenbe@isl . stanford.edu 
boyd@isl . stanford.edu 


Copyright  ©1997  by  Miguel  Sousa  Lobo,  Lieven  Vandenberghe,  and  Stephen  Boyd.  Permis¬ 
sion  to  use,  copy,  modify,  and  distribute  this  Software  for  any  purpose  without  fee  is  hereby 
granted,  provided  that  this  entire  notice  is  included  in  all  copies  of  any  software  which  is 
or  includes  a  copy  or  modification  of  this  software  and  in  all  copies  of  the  supporting  doc¬ 
umentation  for  such  software.  This  software  is  being  provided  “as  is”,  without  any  express 
or  implied  warranty.  In  particular,  the  authors  do  not  make  any  representation  or  warranty 
of  any  kind  concerning  the  merchantability  of  this  software  or  its  fitness  for  any  particular 
purpose. 


1 


1  Introduction 


This  package  contains  software  for  solving  the  second-order  cone  programming  (SOCP)  prob¬ 
lem 

minimize  fTx  .  . 

subject  to  || Apx  +  6* ||  <  cfx  +  di  i  =  1, . . . ,  L. 

The  optimization  variable  is  the  vector  x  G  RTO.  The  problem  data  are  /  G  RTO,  and,  for 

i  =  1 . .  /. .  A, i  G  RNiXTO,  bi  G  R/V\  Cj  G  Rm  and  di  G  R.  The  norm  appearing  in  the 

constraints  is  the  Euclidean  norm,  be.,  ||u||  =  (tTu)1/2.  If  TV*  =  0  we  interpret  the  ?'th 
constraint  as 

0  <  cj x  +  di,  (2) 

i.  e. ,  a  linear  inequality. 

We  refer  to  [LVBL97]  (distributed  with  this  software  package)  for  a  survey  of  applications 
of  second-order  programming. 


1.1  Overview 

The  package  contains: 

•  C-source:  socp .  c  and  socp .  h  contain  a  C-function  socp  that  solves  the  general  second- 
order  cone  programming  problem.  See  §4. 

•  A  Matlab4  interface,  which  allows  the  user  to  call  the  code  from  matlab4.  Compiled 
mex-files  for  some  platforms  are  included;  for  other  platforms  you  can  compile  the 
interface  from  its  source  code  socp_mex.c.  See  §4. 

•  A  Matlab  routine  socp.m.  Although  the  mex-files  can  be  called  directly  from  matlab, 
they  can  also  be  accessed  with  this  routine  that  greatly  simplifies  their  usage.  It 
provides  default  parameters,  diagnostics  and  procedures  for  finding  initial  points.  See 
§2- 

•  Matlab  routine:  socpl  .m.  When  the  problem  has  only  one  constraint,  it  has  an  analytic 
solution,  which  can  be  computed  with  this  function.  See  §3. 

•  Matlab  examples:  we  provide  simple  matlab  routines  that  use  socp.m  to  solve  various 
problems  (some  of  which  are  described  in  the  survey  paper  [LVBL97]).  See  §5. 

•  Documentation,  in  postscript  format:  doc.ps  (this  document). 

•  The  survey  paper  [LVBL97],  in  postscript  format:  socp.ps. 

This  software  package  is  available  at  http :  / /www-isl .  Stanford .  edu/ people/boyd/ group_index .  htm! 
(or  via  anonymous  ftp  from  isl.stanford.edu  in  pub/boyd/). 
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1.2  How  to  use  the  package 

You  can  use  the  package  in  different  ways. 

•  Write  a  C-program  that  calls  the  function  socp  in  socp.c,  and  link  the  code  with 
LAPACK  and  BLAS.  You  should  be  able  to  do  this  on  any  machine  with  a  C  compiler; 
see  §4  for  more  details. 

•  The  easiest  method  (provided  you  have  matlabd)  is  to  use  the  m-file  socp.m.  socp.m 
calls  the  C-function  via  a  mex-file  interface.  Both  socp.m  and  the  mex-files  must  be 
in  your  matlab  path.  See  §2. 

•  If  none  of  the  compiled  mex-files  is  appropriate  for  your  machine,  create  a  mex-file  for 
use  in  matlabd  by  compiling  socp_mex.c  using  cmex,  and  linking  it  with  LAPACK  and 
BLAS  (see  §4  for  more  details). 

1.3  Who  can/should  use  the  code? 

The  code  is  intended  for  researchers  who  want  to  apply  second-order  cone  programming 
(which  includes  quadratically  constrained  quadratic  programming).  Our  aim  is  to  provide  a 
program  that  is  simple,  easy  to  use,  reasonably  efficient,  and  useful  for  small  and  medium¬ 
sized  problems  (say,  up  to  hundreds  of  variables). 

1.4  Caveat 

SOCP  is  a  simple,  general  code,  and  does  not  exploit  problem  structure,  nor  sparsity. 
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2  Matlab  function:  socp.m 

[x , inf o,z,w, hist .time]  =  ... 

socp(f  3A3b3C,d3N3x3z3w3abs_tol 3rel_tol , target 3max_iter 3Nu3out_mode) ; 

Any  of  the  input  parameters  x ,  z,  w,  abs_tol3  rel_tol3  target,  max_iter3  Nu  and 
out_mode,  and/or  the  output  parameters  info3z3w3hist  and  time  can  be  omitted.  They 
must  be  omitted  in  order  from  the  end,  i.e.,  if  a  given  parameter  is  used,  all  the  preceeding 
ones  must  also  be  used.  Alternatively,  any  of  these  parameters  can  be  replaced  by  the  matlab 
empty  list:  [  ] .  In  both  cases,  the  parameter  is  replaced  by  a  default  value  or,  in  the  case  of 
the  initial  points,  an  alternative  procedure  is  used  to  compute  them.  See  detailed  discussion 
of  parameters  below. 

The  shortest  calling  sequence  is: 

x  =  socp(f 3A3b3C,d3N) ; 

2.1  Purpose 

socp  solves  the  second-order  cone  program 

minimize  fTx 

subject  to  \\AiX  +  6* ||  <  c[x  +  d,  i  =  1, . . . ,  L 

and  its  dual 

maximize  —  Y.f= i  (bf  A  +  diwA 
subject  to  £f=i  (Ajzi  +  =  f 

||zj||  <  Wi  i  =  l,...,L 

given  strictly  feasible  primal  and  dual  initial  points. 

2.2  Storage  convention 

The  problem  data  A*  E  RN‘xn,  b,  E  Rv’',  c*  E  Rn,  d,;  E  R,  and  the  dual  variables  z-t,  E  R'¥,:, 
and  Wj,  E  R,  are  stored  as  matrices  and  vectors  A ,  6,  C,  d.  z  and  w,  defined  as  follows: 


'  A1  ' 

'  h  ' 

r  t  i 

c{ 

d\ 

Zl 

W\ 

b  = 

c  = 

d  = 

z  = 

W  = 

A l 

.  bL  _ 

i 

_ i 

d  /, 

.  ZL  _ 

.  U,L  _ 

Note  that  the  socp.m  data  storage  convention  is  different  from  the  convention  used  in  the 
C  and  mex  functions. 


4 


2.3  Input  arguments 

1.  f :  Vector  of  length  n,  specifies  the  primal  objective. 

L 

2.  A:  Matrix  of  size  ^  iV*  by  n,  as  described  above. 

i— 1 

L 

3.  b:  Vector  of  length  ^  A',.  as  described  above. 

2=1 

4.  C:  Matrix  of  size  L  by  n,  as  described  above. 

■5.  d:  Vector  of  length  L,  as  described  above. 

6.  N:  Vector  of  length  L.  The  i-th  element  is  the  number  of  rows  in  .4,. 

7.  x:  Vector  of  length  n.  Primal  starting  point.  Must  be  strictly  feasible.  If  x  is  omitted, 
a  phase  1  method  is  used  to  find  an  initial  primal  point  (See  §s-phasel). 

L 

8.  z:  Vector  of  length  ^  TV*.  Dual  starting  point  (together  with  w). 

2=1 

9.  w:  Vector  of  length  L.  Dual  starting  point  (together  with  z).  z  and  w  must  be  strictly 
dual  feasible.  If  omitted,  a  big -M  constraint  is  added  to  the  SOCP  (See  §s-bigM) 

10.  abs_tol:  Absolute  tolerance.  Default  value:  10-6.  See  discussion  of  convergence 
criteria  below. 

11.  rel_tol:  Relative  tolerance.  Default  value:  10-4.  See  discussion  of  convergence  criteria 
below. 

12.  target:  Target  value,  only  used  if  rel_tol<  0.0.  Default  value:  0.0.  See  discussion  of 
convergence  criteria  below. 

13.  max_iter:  Maximum  number  of  iterations.  Default  value:  100. 

14.  Nu:  The  parameter  u  that  controls  the  rate  of  convergence,  v  >  1.0.  Default  value: 
10. 

15.  out_mode:  Specifies  what  should  be  output  in  hist:  0-  nothing  (default  value):  1-  the 
duality  gap,  for  the  initial  point  and  after  each  iteration;  2-  the  duality  gap  and  the 
deviation  from  centrality,  for  the  initial  point  and  after  each  iteration.  (For  more  on 
duality  gap  and  deviation  from  centrality,  see  [LVBL97]  and  [VB96].) 
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2.4  Output  arguments 

1.  x:  Vector  of  length  n.  The  last  primal  iterate  x. 

2.  info:  String;  possible  values:  ’absolute  accuracy  reached’,  ’relative  accuracy 
reached’,  ’target  value  reached’,  ’target  value  is  unachievable’, 

’maximum  iterations  exceeded’,  ’error’. 

3.  z,  and 

L 

4.  w:  Vectors  of  length  ^  iV*  and  of  length  L.  The  last  dual  iterate  z,w. 

2—1 

■5.  hist:  See  out_mode. 

6.  time:  Vector  with  3  numbers;  performance  statistics:  time(l)  is  user  time  (cpu  time 
used),  time (2)  is  system  time  (time  spent  in  system  calls),  time (3)  is  total  number 
of  iterations  performed. 

2.5  Convergence  criteria 

socp  computes  an  interval  in  which  the  optimal  value  has  been  located.  The  upper  bound  is 
the  primal  objective  p  =  fTx  evaluated  at  the  primal  feasible  solution  stored  in  x;  the  lower 
bound  is  the  dual  objective  cl  =  —  Y%=\  (bf  zi  +  d,ic}j  evaluated  at  the  dual  feasible  solution 
stored  in  z. 

The  length  of  this  interval  {p  —  d)  is  called  the  duality  gap  associated  with  x  and  z 
(see  [LVBL97]). 

The  program  exits  normally  under  five  possible  conditions: 

•  The  maximum  number  of  iterations  is  exceeded. 

•  The  absolute  tolerance  is  reached: 


p  —  d  <  abs_tol. 


•  The  relative  tolerance  is  reached.  The  primal  objective  p  and  the  dual  objective  d  are 
both  positive  and 

p  —  cl 

— - —  <  rel_tol, 

cl 

or  the  primal  and  dual  objective  are  both  negative  and 


p  —  cl 

- <  rel_tol. 

~P 


•  The  target  vcdue  is  reached.  rel_tol  is  negative  and  the  primal  objective  p  is  less  than 

target. 
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•  The  target  value  is  unachievable.  rel_tol  is  negative  and  the  dual  objective  cl  is  greater 

than  target. 

The  target  value  stopping  condition  is  useful  for  feasibility  problems:  socp  stops  when 
it  has  either  found  a  point  with  objective  less  than  the  given  target,  or  has  found  a  dual 
point  with  objective  value  greater  than  target  (which  proves  the  optimal  objective  is  greater 
than  target).  The  target  value  stopping  condition  is  enabled  by  passing  a  negative  relative 
tolerance  argument,  e.g.,  rel_tol  =  —1.0. 


2.6  Caveats 


The  columns  of 


A 

C 


should  be  linearly  independent.  If  this  is  not  the  case,  either 


the  problem  is  not  dual  feasible  (which  means  the  primal  is  unbounded)  or  it  can  be 
reduced  to  another  problem  with  fewer  variables. 


This  reduction  is  done  by  default  in  socp,  but  requires  some  computational  effort.  If, 
due  to  the  nature  of  your  problem,  you  know  that  the  matrix  is  full-rank,  you  can  skip 
the  verification  and  size  reduction  procedure  by  editing  the  appropriate  line  in  the  file 
socp.m  from 


check_rank=l ; 


to 


check_rank=0 ; 

is  such  that  (even  after  the  transformation  above)  its  size  is  close  to  square, 

some  numerical  difficulties  may  arise.  In  this  case  you  should  consider  using  the  dual 
problem  as  the  primal. 


If 


A 

C 


•  The  optimization  procedure  uses  strictly  feasible  primal  and  dual  points,  which  means 
that  the  standard  trick  to  add  equality  constraints  (writing  Ax  =  b  as  the  two  vector 
inequalities  Ax  —  b  >  0,  b  —  Ax  >  0)  will  not  work.  You  have  to  explicitly  eliminate 
equality  constraints. 


2.7  Initial  dual  point:  Big-M 

Note:  If  you  just  want  to  use  socp.m  and  are  not  interested  in  how  it  works,  you  can  skip 
this  and  the  next  sections. 

When  the  dual  variables  z  and  w  are  not  specified  in  the  call  to  socp.m,  a  big-M  procedure 
is  used.  The  original  problem  is  extended  to  include  a  bound  on  the  primal  variable.  The 
dual  of  the  modified  problem  is  such  that  a  strictly  feasible  dual  point  can  always  be  easily 
computed. 

However,  the  solution  to  the  original  problem  may  lie  outside  of  the  added  bound.  In 
this  case,  the  optimization  would  converge  to  a  point  on  the  boundary  of  the  set  defined  by 
the  new  bound,  and  not  to  the  desired  solution.  This  problem  is  circunvented  by  increasing 
the  value  of  the  bound  repeatedly,  to  ensure  it  is  not  active  in  the  solution. 
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The  modified  second-order  cone  program  solved  by  socp.m  is 


minimize  fTx 

subject  to  1 1 A i x  +  b,  \  <  cfx  +  <k  i  =  1, . . . ,  L 
Zil (cfx  +  Ck)  <  M 

and  its  dual  is 

maximize  —  Ef=i  (bf  Zi  +  di(wi  /;))  Mv 
subject  to  E*=i  (Afzi  +  Ci{wi  -«))=/ 

||^i||  <  Wi  i  =  1 . L 

v  >  0. 


A  strictly  feasible  primal  initial  point  is  assumed. 


For 


A 

C 


with  independent  columns,  all  the  primal  variables  Xi  are  bounded  by  the  addi¬ 


tional  constraint.  In  this  case,  the  extra  dual  variables  associated  with  the  added  constraint 
allow  for  the  easy  computation  of  a  strictly  feasible  dual  point. 

M  is  iteratively  increased  to  keep  the  bound  inactive.  Every  BigM_iter  iterations,  if 


L 

M  >  BigM_K  J2(cjx  +  ck). 

i— 1 


the  bound  M  is  changed  to 


L 

M  =  BigM_K  ^{cjx  +  di). 

i— 1 

The  default  values  are  BigM_iter=2  and  BigM_K=2. 

The  solution  returned  by  socp.m  is  to  the  original  problem.  (The  original  convergence 
criteria  are  still  strictly  met,  since  the  duality  gap  of  the  modified  problem  is  an  upper  bound 
on  the  duality  gap  of  the  original  problem.) 

2.8  Initial  primal  point:  Phase  1 

When  the  primal  variable  x  is  not  specified  in  the  call  to  socp.m,  a 
used.  A  feasibility  problem  is  first  solved,  and  its  solution  is  then  used 
the  original  problem. 

The  feasibility  problem  is 

minimize  a 
subject  to  || Apx  +  bj}\  <  cfx  +  cf  +  a  i  =  1, . . L. 

A  strictly  feasible  primal  point  is  always  trivially  found  by  taking  a  large  enough.  As  for 
an  initial  dual  point,  the  big-M  procedure  described  above  is  used.  (In  practice,  socp.m 
constructs  the  modified  problem  and  makes  a  recursive  call  to  itself.) 


phase  1  procedure  is 
as  an  initial  point  for 
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The  convergence  criteria  for  this  problem  are  specified  by  target=0.0  (activated  by 
rel_tol=-l . 0).  socp  will  stop  when  a  <  0;  or  when  it  can  be  shown  that,  subject  to  the 
constraints,  min  a  >  0.  That  is,  it  will  stop  as  soon  as  x  becomes  strictly  feasible  for  the 
original  problem;  or  as  soon  as  it  can  be  shown  that  there  is  no  such  strictly  feasible  point. 

If  a  strictly  feasible  primal  point  is  found,  it  is  then  used  as  an  initial  primal  point  to 
solve  the  original  problem. 

As  a  final  note,  we  must  say  that  these  approaches  ( big-M  and  phase  1)  are  certainly  not 
as  good  as  using  an  infeasible  method.  However  they  have  the  advantage  of  being  simple,  and 
are  effective  enough  for  small  to  medium  sized  problems  (with,  say,  a  few  hundred  variables). 
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3  Matlab  function:  socpl.m 

[x,info ,z,w,fsbl ,bndd]  =  socpl(f ,A)b,c,d) ; 

socpl  computes  an  analytic  solution  to  the  single  constraint  second-order  cone  program 

minimize  fTx 

subject  to  || Ax  +  6||  <  cx  +  d 

and  its  dual 

maximize  —bTz  —  din 
subject  to  ATz  +  cTw  =  / 

||z||  —  w 

Note  that  here  c  is  a  row  vector.  The  problem  is  defined  by  A  E  RTOXn,  b  E  R”\  c  E  Rn 
and  d  E  R.  It’s  variables  are  x  E  Rn,  ,2  G  Rm,  and  w  E  R. 

3.1  Input  arguments 

1.  f :  Vector  of  length  n,  specifies  the  primal  objective. 

2.  A:  Matrix  of  size  rn  by  n. 

3.  b:  Column  vector  of  length  rn. 

4.  c:  Row  vector  of  length  n. 

5.  d:  Scalar. 

3.2  Output  arguments 

1.  x:  Vector  of  length  n.  Primal  solution  x ,  if  one  could  be  computed. 

2.  info:  String;  possible  values:  ’solution  found’,  ’problem  is  infeasible’, 
’problem  is  unbounded’,  ’error’. 

3.  z:  Vector  of  length  rn.  and 

4.  w:  Scalar,  z  and  w  provide  a  dual  solution,  if  one  could  be  computed. 

5.  f  sbl:  Scalar.  1  if  the  problem  is  feasible,  0  otherwise. 

6.  bndd:  Scalar.  1  if  the  problem  has  a  bounded  solution,  0  otherwise. 
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4  C-function:  socp.c 

The  main  routine  is  a  C-function  in  socp.c: 


int  socp(  int  L,  int  *N,  int  n,  double  *f3  double  *A3  double  *b3 
double  *x,  double  *z3 

double  abs_tol3  double  rel_tol3  double  target,  int  *iter, 
double  Nu3 

int  *info,  int  out_mode,  double  *hist3 
double  *dblwork,  int  *intwork  ) 

To  compute  the  amount  of  workspace  required  (dblwork,  ptrwork  and  intwork).  another 
C-function  is  provided  in  socp.c.  This  should  be  called  before  socp: 

void  socp_getwork(  int  L,  int  *N3  int  n,  int  max_iter3  int  out_mode3 

int  *mhist3  int  *nhist3 
int  *ndbl ,  int  *nint  ) 


4.1  Purpose 

socp  solves  the  second-order  cone  program 

minimize  fTx 

subject  to  ||  AiX  +  b,  \  <  cjx  +  <A,  i  =  1, . . . ,  L 

and  its  dual 

maximize  —  Y%=\  (bj z%  +  ckwi) 
subject  to  Ef=  i  (  i/  I  +  CjtOi)  =  / 

||^i||  <  Wi  i  =  l,...,L 

given  strictly  feasible  primal  and  dual  initial  points. 


4.2  Storage  convention 

The  sizes  of  the  constraints  are  given  in  a  vector  N,  of  length  L.  iV,;  defines  the  dimensionality 
of  the  /-tli  conic  constraint  ,  i.e.,  the  dimension  of  the  argument  of  the  norm  plus  one.  Hence, 
,1,  €  RiV‘  n  :"',  b,  €  R;  V;  p;.  a  e  Rn,  d,  G  R,  Zi  e  RiV;  P;.  and  mt  €  R. 

For  the  C  and  mex  functions,  A,  6,  and  the  dual  variable  z,  are  conventioned  to  be 


A\ 

'  h  ' 

ft 

T 

C1 

rl| 

Wi 

A  = 

al 

b  = 

bL 

z  = 

Zl 

i 

_ 1 

<h- 

.  WL  . 

Note  that  this  is  different  from  the  data  storage  convention  used  in  socp.m. 
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4.3  Arguments  for  socp 

1.  L  (integer,  in):  The  number  of  conic  constraints. 

2.  N  (pointer  to  integer,  in):  Array  of  length  L.  Thez-th  element  is  the  dimension  of  z-tli 
constraint. 

3.  n  (integer,  in):  The  number  of  variables. 

4.  f  (pointer  to  double,  in):  Array  of  length  n,  specifies  the  primal  objective. 

5.  A  (pointer  to  double,  in):  Array  of  length  ran.  as  described  above  (storage  is  one 
column  at  a  time,  i.e .,  element  t.  j  of  matrix  is  entry  i  +  rn  *  (j  —  1)  of  array). 

6.  b  (pointer  to  double,  in):  Array  of  length  m,  as  described  above. 

7.  x  (pointer  to  double,  in/out):  Array  of  length  n.  On  entry,  a  strictly  feasible  primal 
point.  On  exit,  last  primal  iterate  x.  Meaning  depends  on  stopping  criterion  (solution 
to  primal  problem,  if  appropriate). 

8.  z  (pointer  to  double,  in/out):  Array  of  length  m.  On  entry,  a  strictly  feasible  dual 
point,  as  described  above.  On  exit,  last  dual  iterate  z.  Meaning  depends  on  stopping 
criterion  (solution  to  dual  problem,  if  appropriate). 

9.  abs_tol  (double,  in):  Absolute  tolerance.  See  discussion  of  convergence  criteria 
below. 

10.  rel_tol  (double,  in):  Relative  tolerance.  See  discussion  of  convergence  criteria  below. 

11.  target  (double,  in):  Target  value,  only  used  if  rel_tol<  0.0.  See  discussion  of 
convergence  criteria  below. 

12.  iter  (pointer  to  integer,  in/out):  On  entry,  *iter  is  the  maximum  number  of 
iterations,  socp  is  aborted  if  more  are  required  for  convergence.  On  exit,  *iter  is  the 
number  of  iterations  performed. 

13.  Nu  (double,  in):  The  parameter  v  that  controls  the  rate  of  convergence,  v  >  1.0. 
Recommended  range:  5.0-50.0. 

14.  info  (pointer  to  integer,  out):  Stopping  criteria  that  caused  exiting;  0-  error;  1- 
absolute  tolerance;  2-  relative  tolerance;  3-  target  value  (was  achieved  or  is 
unachievable);  4-maximum  iterations. 

15.  out_mode  (integer,  in):  Specifies  what  should  be  output  in  hist;  0-  nothing;  1- 
duality  gap,  for  initial  point  and  after  each  iteration;  2-  duality  gap  and  deviation 
from  centrality,  for  initial  point  and  after  each  iteration. 

16.  hist  (pointer  to  double,  out):  Array  of  length  greater  or  equal  to  the  product  of 
*mhist  by  *nhist.  See  out_mode. 
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17.  dblwork  (pointer  to  double):  Array  of  doubles  for  workspace.  See  socp_getwork  for 
length. 

18.  intwork  (pointer  to  integer):  Array  of  integers  for  workspace.  See  socp_getwork  for 
length. 

socp  returns  0  if  it  exited  due  to  any  of  the  stopping  criteria,  or  an  error  code  from  the 
LAPACK  routine  dgelss. 

4.4  Arguments  for  socp_getwork 

1.  n  (integer,  in):  Use  same  value  as  in  socp. 

2.  L  (integer,  in):  Use  same  value  as  in  socp. 

3.  N  (pointer  to  integer,  in):  Use  same  values  as  in  socp. 

4.  max_iter  (integer,  in):  Use  entry  value  of  *iter  in  socp. 

5.  out_mode  (integer,  in):  Use  same  value  as  in  socp. 

6.  mhist  (pointer  to  integer,  out),  and 

7.  nhist  (pointer  to  integer,  out):  The  required  length  of  hist  in  number  of  doubles  is 

the  product  of  mhist  by  nhist.  (These  are  included  to  make  it  easy  to  edit  the  code  to 
send  other  output  through  hist;  currently  the  values  are  always  1  and  max_iter  +  1, 

since  hist  returns  a  vector  with  the  duality  gap  for  the  initial  point  and  after  each 

iteration.) 

8.  ndbl  (pointer  to  integer,  out):  Returns  required  length  for  dblwork,  in  number  of 
doubles. 

9.  nint  (pointer  to  integer,  out):  Returns  required  length  for  intwork,  in  number  of 
integers. 

4.5  Convergence  criteria 

(This  is  the  same  as  for  socp.mj 

socp  computes  an  interval  in  which  the  optimal  value  has  been  located.  The  upper 
bound  is  the  primal  objective  p  =  fTx  evaluated  at  the  primal  feasible  solution  stored  in  x; 
the  lower  bound  is  the  dual  objective  d  =  —  Yn=\  (bf A  +  diWyj  evaluated  at  the  dual  feasible 
solution  stored  in  z. 

The  length  of  this  interval  (p  —  cl)  is  called  the  duality  gap  associated  with  x  and  z 
(see  [LVBL97]). 

The  program  exits  normally  under  five  possible  conditions: 

•  The  maximum  number  of  iterations  is  exceeded. 
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The  absolute  tolerance  is  reached: 


p  —  cl  <  abs_tol. 


•  The  relative  tolerance  is  reached.  The  primal  objective  p  and  the  dual  objective  d  are 
both  positive  and 

p  —  cl 

— - —  <  rel_tol, 

cl 

or  the  primal  and  dual  objective  are  both  negative  and 


p  —  cl 

- <  rel_tol. 

~P 


•  The  target  value  is  reached.  rel_tol  is  negative  and  the  primal  objective  p  is  less  than 

target. 

•  The  target  vcdue  is  unachievable.  rel_tol  is  negative  and  the  dual  objective  cl  is  greater 

than  target. 

The  target  value  stopping  condition  is  useful  for  feasibility  problems:  socp  stops  when 
it  has  either  found  a  point  with  objective  less  than  the  given  target,  or  has  found  a  dual 
point  with  objective  value  greater  than  target  (which  proves  the  optimal  objective  is  greater 
than  target).  The  target  value  stopping  condition  is  enabled  by  passing  a  negative  relative 
tolerance  argument,  e.g.,  rel_tol  =  —1.0. 


4.6  Caveats 

•  Little  or  no  effort  in  parameter  validation  is  made  in  these  functions,  so  some  surprises 
are  to  be  expected  if  they  are  used  without  care.  The  mex  and  m-files  include  parameter 
validation. 

•  For  socp.c,  the  columns  of  A  must  be  linearly  independent.  If  this  is  not  the  case, 
either  the  problem  is  not  dual  feasible  (which  means  the  primal  is  unbounded)  or  it 
can  be  reduced  to  one  with  fewer  variables.  This  reduction  is  done  automatically  in 
the  m-files. 

•  If  A  is  such  that  its  size  is  close  to  square,  some  numerical  difficulties  may  arise.  In 
this  case  you  should  consider  using  the  dual  problem  as  the  primal. 

•  The  big-M  and  phase  1  procedures  are  not  implemented  in  C;  a  strictly  feasible  primal 
and  dual  initial  point  must  be  provided. 

•  The  optimization  procedure  uses  strictly  feasible  primal  and  dual  points,  which  means 
that  standard  trick  of  adding  equality  constraints  (writing  Ax  =  b  as  the  two  vector 
inequalities  Ax  —  b  >  0,  b  —  Ax  >  0)  will  not  work.  You  have  to  explicitly  eliminate 
equality  constraints. 
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4.7  Compiling 

The  source  code  for  the  function  socp  is  in  the  files  socp.c  and  socp.h.  It  is  written  in 
ansi-C  with  calls  to  LAPACK  and  BLAS. 

LAPACK  can  be  obtained  via  Netlib  (http://www.netlib.org  or  anonymous  ftp  at 
ftp.netlib.org).  A  set  of  optimized  BLAS-routines  should  be  supplied  by  your  computer 
vendor.  A  non-optimized  version  can  also  be  obtained  from  Netlib. 

You  can  also  create  a  mex-file  for  use  in  matlabl  by  compiling  socp_mex.c  using  cmex, 
and  linking  it  with  LAPACK  and  BLAS.  The  details  of  compiling  the  code  vary  with  the 
compiler  used  and  other  factors  such  as  where  various  libraries  are  located.  We  have  provided 
a  sample  Makefile  in  the  source  directory.  In  order  to  compile  the  mex-interface,  edit  this 
Makefile  as  indicated,  and  type  make. 

Executable  mex-files  are  provided  for  some  platforms. 
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5  Example  mat  lab  files 

Most  of  the  examples  provided  can  be  easily  be  understood  by  reading  the  initial  comments 
in  the  files.  We  provide  here  some  additional  comments  on  some  of  the  examples. 


5.1  grasp 

The  files  grasp.l.m  and  grasp_robust .m  in  the  directory  examples/grasp  demonstrate 
the  grasping  problem  described  in  [LVBL97]  (this  paper  includes  a  diagram  of  the  problem 
geometry).  The  solution  is  computed  by  calling  the  script  grasp_polyhedron.m.  This  is  a 
general  script  that  can  be  used  to  treat  other  grasping  problems. 


•  grasp_polyhedron.m 

grasp_polyhedron.m  is  a  script  that  transforms  a  robust  grasping  problem  into  an 
SOCP.  The  objective  is:  find  the  maximum  scaling  factor  K  such  that  there  is  a  grasp 
that  stabilizes  all  the  force  and  torque  pairs 


^ext  —  ft  Tj  +  Fb 
T^t  =  KTi  +  Th 


i  =  1 


Such  a  grasp  will  be  stable  for  any  force  and  torque  in  the  convex  hull  of  the  Fext  and 
Text  a  polyhedron  in  R6). 

grasp_polyhedron.m  uses  as  input  the  following  variables,  that  must  be  pre-defined 
(n  is  the  number  of  contact  points,  and  m  is  the  number  of  vertices  in  the  force/torque 
polyhedron): 

1.  p  (matrix,  3  by  n):  position  relative  to  center  of  mass  for  each  contact  point. 

2.  u  (matrix,  3  by  n):  direction  of  force  for  each  contact  point. 

3.  v  (matrix,  3  by  n):  inward  normal  to  the  surface  for  each  contact  point. 

4.  Fa  (matrix  3  by  m):  vertices  of  polyhedron,  force  component  (relative  to  its  center 

Fb). 

5.  Ta  (matrix,  3  by  m):  vertices  of  polyhedron,  torque  component  (relative  to  its 
center  Tb). 

6.  fmax  (scalar):  upper  bound  on  forces  f. 

The  script  writes  its  output  in  the  following  variables: 

1.  K  (scalar):  the  maximum  scaling  factor  for  which  there  is  a  stable  grasp. 

2.  f  (vector,  length  n):  the  grasping  forces  in  the  corresponding  to  the  max- it  grasp, 
for  each  contact  point. 

3.  F  (matrix,  3  by  nm ):  the  friction  forces  corresponding  to  the  max- it  grasp,  for 
each  contact  point,  and  for  each  external  loading. 
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Note  that  the  script  uses  temporary  variables,  which  may  overwrite  variables  in  the 
workspace  with  the  same  name. 

•  grasp_l.m 

grasp_l  .m  is  a  script  that  defines  the  simple  example  described  in  [LVBL97],  and  calls 
grasp_polyhedron.m.  The  problem  is  solved  for  different  friction  coefficients  //.  For 
each  value  of  //  we  compute  the  maximum  torque  for  which  there  is  a  stable  grasp, 
with  upper  bounds  on  the  grasping  forces.  It  produces  the  figure  shown  in  the  paper. 

•  grasp_robust .m 

grasp_robust .  m  is  a  script  that  demonstrates  robust  grasping,  with  the  approximation 
of  a  ball  of  possible  forces  by  a  polyhedron  with  m  vertices.  As  m  increases,  the  poly¬ 
hedron  converges  to  the  ball  (with  probability  one),  and  the  results  of  the  optimization 
also  converge. 

5.2  matrix-frac 

The  script  mf  .m,  in  the  directory  examples/matrix-frac,  demonstrates  the  matrix-fractional 
programming  problem  described  in  [LVBL97]. 

A  particularity  of  this  problem  is  that  the  solution  to  the  original  problem  is  obtained 
from  the  dual  variable: 

xp=zeros (p, 1) ; 
for  i=l:p, 

xp(i)=l/2*(z((i+l)* (n+1) ) -w(i+l) ) ; 
end; 

The  script  repeats  random  problems  over  a  range  of  sizes,  and  several  times  for  each  size 
to  obtain  statistical  information  on  the  performance  of  the  optimization  method. 

5.3  lp 

The  script  lp.m,  in  examples/lp,  generates  and  solves  random  linear  programs.  It  has  the 
following  possibilities: 

•  Test  a  range  of  Mu;  this  can  be  controled  by  editing  the  line 
Nu= [10  20  100]; 

•  Test  a  range  of  problem  sizes  n;  defined  by  the  line 
for  n=20:20:100, 

•  For  each  value  of  Mu  and  n,  repeat  to  obtain  statistical  measures,  average  and  variance; 
this  is  defined  by  the  line 

for  k=l :  11 , . 
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